Carga y limpieza preliminar de los datos

Los datos que se van a analizar proceden de Kaggle - Version 91 - 29 de abril de 2020

insertando un chunk de Python

#{python} #import pandas as pd #datos = pd.read_csv("covid_19_clean_complete.csv") #datos.head(10) #variable_chunk_python = 2 #

insertando un chunk de R y usando ‘reticulate’ que hace de conexion entre R y Python

#pd <- import("pandas")
#datos <- pd$read_csv("covid_19_clean_complete.csv")
#kable(head(datos, 10))
# Como llamar a una variable creada en un chunk de python
# variable_chunk_r = variable_chunk_python * 3 #--> esto da error
#variable_chunk_r = py$variable_chunk_python * 3 #--> esto es Ok, poniendo 'py$' antes de la variable

importando los datos nativamente en R

datos <-read.csv("covid_19_clean_complete.csv", stringsAsFactors = FALSE)
#kable(head(datos, 10))
# utilizando la libreria 'tidyverse'
datos %>% head(10) %>% kable()
Province.State Country.Region Lat Long Date Confirmed Deaths Recovered
Afghanistan 33.0000 65.0000 1/22/20 0 0 0
Albania 41.1533 20.1683 1/22/20 0 0 0
Algeria 28.0339 1.6596 1/22/20 0 0 0
Andorra 42.5063 1.5218 1/22/20 0 0 0
Angola -11.2027 17.8739 1/22/20 0 0 0
Antigua and Barbuda 17.0608 -61.7964 1/22/20 0 0 0
Argentina -38.4161 -63.6167 1/22/20 0 0 0
Armenia 40.0691 45.0382 1/22/20 0 0 0
Australian Capital Territory Australia -35.4735 149.0124 1/22/20 0 0 0
New South Wales Australia -33.8688 151.2093 1/22/20 0 0 0

Estructura de los datos

str(datos)
## 'data.frame':    25676 obs. of  8 variables:
##  $ Province.State: chr  "" "" "" "" ...
##  $ Country.Region: chr  "Afghanistan" "Albania" "Algeria" "Andorra" ...
##  $ Lat           : num  33 41.2 28 42.5 -11.2 ...
##  $ Long          : num  65 20.17 1.66 1.52 17.87 ...
##  $ Date          : chr  "1/22/20" "1/22/20" "1/22/20" "1/22/20" ...
##  $ Confirmed     : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Deaths        : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Recovered     : int  0 0 0 0 0 0 0 0 0 0 ...
colnames(datos) = c("Provincia_Estado", # Cualitativo
                    "Pais_Region", # Cualitativo
                    "Latitud", # Norte+ o Sur- Cuantitativo
                    "Longitud", # Este+ u Oeste- Cuantitativo
                    "Fecha", # Ordinal
                    "Casos_Confirmados", # Cuantitativo
                    "Casos_Muertos", # Cuantitativo
                    "Casos_Recuperados" # Cuantitativo
                    )
datos %>% head() %>% kable() #%>% kable_styling()
Provincia_Estado Pais_Region Latitud Longitud Fecha Casos_Confirmados Casos_Muertos Casos_Recuperados
Afghanistan 33.0000 65.0000 1/22/20 0 0 0
Albania 41.1533 20.1683 1/22/20 0 0 0
Algeria 28.0339 1.6596 1/22/20 0 0 0
Andorra 42.5063 1.5218 1/22/20 0 0 0
Angola -11.2027 17.8739 1/22/20 0 0 0
Antigua and Barbuda 17.0608 -61.7964 1/22/20 0 0 0
  • Cualitativas se convierten con ‘factor’ o bien ‘as.factor’.
  • Ordinales se convierten con ‘ordered’.
  • Cuantitativos se convierten con ‘as.numeric’ o ‘as.integer’.

Cuando se usa en read.csv hay un parámetros llamado ‘stringsAsFactors’ que por defecto es True, pero si lo ponemos a False, no convertirá las columnas string en factores, se quedaran como ‘chr’. Por lo tanto luego sería necesario convertirlas con factor(), las fechas tambien se importan como ‘chr’.

# importado con stringsAsFactors = False
#datos$Provincia_Estado = factor(datos$Provincia_Estado)
#datos$Pais_Region = factor(datos$Pais_Region)
# si importamos magrittr podemos usar %<>% y cambiar estas dos últimas instrucciones
datos$Provincia_Estado %<>% factor() # Primero pasa el valor a la función (dcha) y luego se devulve convertido al propio valor (izda)
datos$Pais_Region %<>% factor()
# Al cargar los datos con stringsAsFactors = False, la fecha se importa como 'chr' y hay que convertirla indicando el formato 'm/d/y' en que vienen los datos en el fichero
#datos$Fecha %<>% as.Date(format="%m/%d/&y")
# Importamos la libreria 'lubridate' que sobreescribe cosas del paquete estandar, y entre ellas es el objeto date y ya no funciona 'as.Date'
datos$Fecha %<>% mdy()
str(datos)
## 'data.frame':    25676 obs. of  8 variables:
##  $ Provincia_Estado : Factor w/ 81 levels "","Alberta","Anguilla",..: 1 1 1 1 1 1 1 1 6 49 ...
##  $ Pais_Region      : Factor w/ 185 levels "Afghanistan",..: 1 2 3 4 5 6 7 8 9 9 ...
##  $ Latitud          : num  33 41.2 28 42.5 -11.2 ...
##  $ Longitud         : num  65 20.17 1.66 1.52 17.87 ...
##  $ Fecha            : Date, format: "2020-01-22" "2020-01-22" ...
##  $ Casos_Confirmados: int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Casos_Muertos    : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Casos_Recuperados: int  0 0 0 0 0 0 0 0 0 0 ...

\[Confirmados = Muertos + Recuperados + Enfermos \]

# Añadimos una nueva columna calculada según la formula anterior, dado que no tenemos los Enfermos
# Con la sintaxis %<>% los datos van a la dcha se añade la columna y vuelven a la izda
datos %<>%
  mutate(Casos_Enfermos = Casos_Confirmados - Casos_Muertos - Casos_Recuperados)
datos %>% 
  filter(Casos_Confirmados > 10000) %>% 
  head(10) %>% 
  kable()
Provincia_Estado Pais_Region Latitud Longitud Fecha Casos_Confirmados Casos_Muertos Casos_Recuperados Casos_Enfermos
Hubei China 30.9756 112.2707 2020-02-02 11177 350 295 10532
Hubei China 30.9756 112.2707 2020-02-03 13522 414 386 12722
Hubei China 30.9756 112.2707 2020-02-04 16678 479 522 15677
Hubei China 30.9756 112.2707 2020-02-05 19665 549 633 18483
Hubei China 30.9756 112.2707 2020-02-06 22112 618 817 20677
Hubei China 30.9756 112.2707 2020-02-07 24953 699 1115 23139
Hubei China 30.9756 112.2707 2020-02-08 27100 780 1439 24881
Hubei China 30.9756 112.2707 2020-02-09 29631 871 1795 26965
Hubei China 30.9756 112.2707 2020-02-10 31728 974 2222 28532
Hubei China 30.9756 112.2707 2020-02-11 33366 1068 2639 29659 
# Comprobar si hay errores en los datos y no concuerdan, sales Casos_Enfermos negativos
datos %>%
  filter(Casos_Enfermos < 0) %>%
  arrange(Provincia_Estado, Fecha) %>%
  kable()
Provincia_Estado Pais_Region Latitud Longitud Fecha Casos_Confirmados Casos_Muertos Casos_Recuperados Casos_Enfermos
Diamond Princess Canada 0.0000 0.0000 2020-03-22 0 1 0 -1
Diamond Princess Canada 0.0000 0.0000 2020-03-23 0 1 0 -1
Diamond Princess Canada 0.0000 0.0000 2020-03-24 0 1 0 -1
Diamond Princess Canada 0.0000 0.0000 2020-03-25 0 1 0 -1
Diamond Princess Canada 0.0000 0.0000 2020-03-26 0 1 0 -1
Diamond Princess Canada 0.0000 0.0000 2020-03-27 0 1 0 -1
Diamond Princess Canada 0.0000 0.0000 2020-03-28 0 1 0 -1
Diamond Princess Canada 0.0000 0.0000 2020-03-29 0 1 0 -1
Diamond Princess Canada 0.0000 0.0000 2020-03-30 0 1 0 -1
Diamond Princess Canada 0.0000 0.0000 2020-03-31 0 1 0 -1
Diamond Princess Canada 0.0000 0.0000 2020-04-01 0 1 0 -1
Diamond Princess Canada 0.0000 0.0000 2020-04-02 0 1 0 -1
Diamond Princess Canada 0.0000 0.0000 2020-04-03 0 1 0 -1
Diamond Princess Canada 0.0000 0.0000 2020-04-04 0 1 0 -1
Diamond Princess Canada 0.0000 0.0000 2020-04-05 0 1 0 -1
Diamond Princess Canada 0.0000 0.0000 2020-04-06 0 1 0 -1
Diamond Princess Canada 0.0000 0.0000 2020-04-07 0 1 0 -1
Diamond Princess Canada 0.0000 0.0000 2020-04-08 0 1 0 -1
Diamond Princess Canada 0.0000 0.0000 2020-04-09 0 1 0 -1
Diamond Princess Canada 0.0000 0.0000 2020-04-10 0 1 0 -1
Diamond Princess Canada 0.0000 0.0000 2020-04-11 0 1 0 -1
Diamond Princess Canada 0.0000 0.0000 2020-04-12 0 1 0 -1
Diamond Princess Canada 0.0000 0.0000 2020-04-13 0 1 0 -1
Diamond Princess Canada 0.0000 0.0000 2020-04-14 0 1 0 -1
Diamond Princess Canada 0.0000 0.0000 2020-04-15 0 1 0 -1
Diamond Princess Canada 0.0000 0.0000 2020-04-16 0 1 0 -1
Diamond Princess Canada 0.0000 0.0000 2020-04-17 0 1 0 -1
Diamond Princess Canada 0.0000 0.0000 2020-04-18 0 1 0 -1
Diamond Princess Canada 0.0000 0.0000 2020-04-19 0 1 0 -1
Diamond Princess Canada 0.0000 0.0000 2020-04-20 0 1 0 -1
Diamond Princess Canada 0.0000 0.0000 2020-04-21 0 1 0 -1
Diamond Princess Canada 0.0000 0.0000 2020-04-22 0 1 0 -1
Diamond Princess Canada 0.0000 0.0000 2020-04-23 0 1 0 -1
Diamond Princess Canada 0.0000 0.0000 2020-04-24 0 1 0 -1
Diamond Princess Canada 0.0000 0.0000 2020-04-25 0 1 0 -1
Diamond Princess Canada 0.0000 0.0000 2020-04-26 0 1 0 -1
Diamond Princess Canada 0.0000 0.0000 2020-04-27 0 1 0 -1
Diamond Princess Canada 0.0000 0.0000 2020-04-28 0 1 0 -1
Hainan China 19.1959 109.7453 2020-03-24 168 6 168 -6
Hainan China 19.1959 109.7453 2020-03-25 168 6 168 -6
Hainan China 19.1959 109.7453 2020-03-26 168 6 168 -6
Hainan China 19.1959 109.7453 2020-03-27 168 6 168 -6
Hainan China 19.1959 109.7453 2020-03-28 168 6 168 -6
Hainan China 19.1959 109.7453 2020-03-29 168 6 168 -6
Hainan China 19.1959 109.7453 2020-03-30 168 6 168 -6
Hainan China 19.1959 109.7453 2020-03-31 168 6 168 -6
Hainan China 19.1959 109.7453 2020-04-01 168 6 168 -6 
# Chequear que pasa con los datos erroneos de la provincia 'Hainan'
datos %>%
  filter(Provincia_Estado == "Hainan") %>%
  kable()
Provincia_Estado Pais_Region Latitud Longitud Fecha Casos_Confirmados Casos_Muertos Casos_Recuperados Casos_Enfermos
Hainan China 19.1959 109.7453 2020-01-22 4 0 0 4
Hainan China 19.1959 109.7453 2020-01-23 5 0 0 5
Hainan China 19.1959 109.7453 2020-01-24 8 0 0 8
Hainan China 19.1959 109.7453 2020-01-25 19 0 0 19
Hainan China 19.1959 109.7453 2020-01-26 22 0 0 22
Hainan China 19.1959 109.7453 2020-01-27 33 1 0 32
Hainan China 19.1959 109.7453 2020-01-28 40 1 0 39
Hainan China 19.1959 109.7453 2020-01-29 43 1 0 42
Hainan China 19.1959 109.7453 2020-01-30 46 1 1 44
Hainan China 19.1959 109.7453 2020-01-31 52 1 1 50
Hainan China 19.1959 109.7453 2020-02-01 62 1 1 60
Hainan China 19.1959 109.7453 2020-02-02 64 1 4 59
Hainan China 19.1959 109.7453 2020-02-03 72 1 4 67
Hainan China 19.1959 109.7453 2020-02-04 80 1 5 74
Hainan China 19.1959 109.7453 2020-02-05 99 1 5 93
Hainan China 19.1959 109.7453 2020-02-06 106 1 8 97
Hainan China 19.1959 109.7453 2020-02-07 117 2 10 105
Hainan China 19.1959 109.7453 2020-02-08 124 2 14 108
Hainan China 19.1959 109.7453 2020-02-09 131 3 19 109
Hainan China 19.1959 109.7453 2020-02-10 138 3 19 116
Hainan China 19.1959 109.7453 2020-02-11 144 3 20 121
Hainan China 19.1959 109.7453 2020-02-12 157 4 27 126
Hainan China 19.1959 109.7453 2020-02-13 157 4 30 123
Hainan China 19.1959 109.7453 2020-02-14 159 4 43 112
Hainan China 19.1959 109.7453 2020-02-15 162 4 39 119
Hainan China 19.1959 109.7453 2020-02-16 162 4 52 106
Hainan China 19.1959 109.7453 2020-02-17 163 4 59 100
Hainan China 19.1959 109.7453 2020-02-18 163 4 79 80
Hainan China 19.1959 109.7453 2020-02-19 168 4 84 80
Hainan China 19.1959 109.7453 2020-02-20 168 4 86 78
Hainan China 19.1959 109.7453 2020-02-21 168 4 95 69
Hainan China 19.1959 109.7453 2020-02-22 168 4 104 60
Hainan China 19.1959 109.7453 2020-02-23 168 5 106 57
Hainan China 19.1959 109.7453 2020-02-24 168 5 116 47
Hainan China 19.1959 109.7453 2020-02-25 168 5 124 39
Hainan China 19.1959 109.7453 2020-02-26 168 5 129 34
Hainan China 19.1959 109.7453 2020-02-27 168 5 131 32
Hainan China 19.1959 109.7453 2020-02-28 168 5 133 30
Hainan China 19.1959 109.7453 2020-02-29 168 5 148 15
Hainan China 19.1959 109.7453 2020-03-01 168 5 149 14
Hainan China 19.1959 109.7453 2020-03-02 168 5 151 12
Hainan China 19.1959 109.7453 2020-03-03 168 5 155 8
Hainan China 19.1959 109.7453 2020-03-04 168 5 158 5
Hainan China 19.1959 109.7453 2020-03-05 168 6 158 4
Hainan China 19.1959 109.7453 2020-03-06 168 6 158 4
Hainan China 19.1959 109.7453 2020-03-07 168 6 158 4
Hainan China 19.1959 109.7453 2020-03-08 168 6 159 3
Hainan China 19.1959 109.7453 2020-03-09 168 6 159 3
Hainan China 19.1959 109.7453 2020-03-10 168 6 159 3
Hainan China 19.1959 109.7453 2020-03-11 168 6 159 3
Hainan China 19.1959 109.7453 2020-03-12 168 6 160 2
Hainan China 19.1959 109.7453 2020-03-13 168 6 160 2
Hainan China 19.1959 109.7453 2020-03-14 168 6 160 2
Hainan China 19.1959 109.7453 2020-03-15 168 6 160 2
Hainan China 19.1959 109.7453 2020-03-16 168 6 161 1
Hainan China 19.1959 109.7453 2020-03-17 168 6 161 1
Hainan China 19.1959 109.7453 2020-03-18 168 6 161 1
Hainan China 19.1959 109.7453 2020-03-19 168 6 161 1
Hainan China 19.1959 109.7453 2020-03-20 168 6 161 1
Hainan China 19.1959 109.7453 2020-03-21 168 6 161 1
Hainan China 19.1959 109.7453 2020-03-22 168 6 161 1
Hainan China 19.1959 109.7453 2020-03-23 168 6 161 1
Hainan China 19.1959 109.7453 2020-03-24 168 6 168 -6
Hainan China 19.1959 109.7453 2020-03-25 168 6 168 -6
Hainan China 19.1959 109.7453 2020-03-26 168 6 168 -6
Hainan China 19.1959 109.7453 2020-03-27 168 6 168 -6
Hainan China 19.1959 109.7453 2020-03-28 168 6 168 -6
Hainan China 19.1959 109.7453 2020-03-29 168 6 168 -6
Hainan China 19.1959 109.7453 2020-03-30 168 6 168 -6
Hainan China 19.1959 109.7453 2020-03-31 168 6 168 -6
Hainan China 19.1959 109.7453 2020-04-01 168 6 168 -6
Hainan China 19.1959 109.7453 2020-04-02 168 6 162 0
Hainan China 19.1959 109.7453 2020-04-03 168 6 162 0
Hainan China 19.1959 109.7453 2020-04-04 168 6 162 0
Hainan China 19.1959 109.7453 2020-04-05 168 6 162 0
Hainan China 19.1959 109.7453 2020-04-06 168 6 162 0
Hainan China 19.1959 109.7453 2020-04-07 168 6 162 0
Hainan China 19.1959 109.7453 2020-04-08 168 6 162 0
Hainan China 19.1959 109.7453 2020-04-09 168 6 162 0
Hainan China 19.1959 109.7453 2020-04-10 168 6 162 0
Hainan China 19.1959 109.7453 2020-04-11 168 6 162 0
Hainan China 19.1959 109.7453 2020-04-12 168 6 162 0
Hainan China 19.1959 109.7453 2020-04-13 168 6 162 0
Hainan China 19.1959 109.7453 2020-04-14 168 6 162 0
Hainan China 19.1959 109.7453 2020-04-15 168 6 162 0
Hainan China 19.1959 109.7453 2020-04-16 168 6 162 0
Hainan China 19.1959 109.7453 2020-04-17 168 6 162 0
Hainan China 19.1959 109.7453 2020-04-18 168 6 162 0
Hainan China 19.1959 109.7453 2020-04-19 168 6 162 0
Hainan China 19.1959 109.7453 2020-04-20 168 6 162 0
Hainan China 19.1959 109.7453 2020-04-21 168 6 162 0
Hainan China 19.1959 109.7453 2020-04-22 168 6 162 0
Hainan China 19.1959 109.7453 2020-04-23 168 6 162 0
Hainan China 19.1959 109.7453 2020-04-24 168 6 162 0
Hainan China 19.1959 109.7453 2020-04-25 168 6 162 0
Hainan China 19.1959 109.7453 2020-04-26 168 6 162 0
Hainan China 19.1959 109.7453 2020-04-27 168 6 162 0
Hainan China 19.1959 109.7453 2020-04-28 168 6 162 0 
# Hay un error entre una fechas, porque luego se corrige
datos %>%
  filter(Provincia_Estado == "Hainan", Casos_Enfermos < 0) %>%
  mutate(Casos_Recuperados = Casos_Recuperados + Casos_Enfermos,
         Casos_Enfermos = 0) %>%
  kable()
Provincia_Estado Pais_Region Latitud Longitud Fecha Casos_Confirmados Casos_Muertos Casos_Recuperados Casos_Enfermos
Hainan China 19.1959 109.7453 2020-03-24 168 6 162 0
Hainan China 19.1959 109.7453 2020-03-25 168 6 162 0
Hainan China 19.1959 109.7453 2020-03-26 168 6 162 0
Hainan China 19.1959 109.7453 2020-03-27 168 6 162 0
Hainan China 19.1959 109.7453 2020-03-28 168 6 162 0
Hainan China 19.1959 109.7453 2020-03-29 168 6 162 0
Hainan China 19.1959 109.7453 2020-03-30 168 6 162 0
Hainan China 19.1959 109.7453 2020-03-31 168 6 162 0
Hainan China 19.1959 109.7453 2020-04-01 168 6 162 0

Análisis geográfico

# Filtrado directo rectangulo aprox. de Europa
#datos_europa = datos[datos$Latitud > 38 & datos$Longitud > -25 & datos$Longitud < 30 , ]
#datos_europa <- datos[datos$Latitud > 38 & datos$Longitud > -25 & datos$Longitud < 30 , ]
# Filtrao de otra manera con la función filter
datos_europa = datos %>%
  filter(Latitud > 38, between(Longitud, -25, 30))

nrow(datos_europa)
## [1] 4410
table(datos_europa$Pais_Region) %>% nrow()
## [1] 185
# ejecutando (Knit) aparecen 185 paises, muchos con frecuencia 0, y que no están en datos_europa ???
frec_europa = as.data.frame(table(datos_europa$Pais_Region))
frec_europa %>% kable()
Var1 Freq
Afghanistan 0
Albania 98
Algeria 0
Andorra 98
Angola 0
Antigua and Barbuda 0
Argentina 0
Armenia 0
Australia 0
Austria 98
Azerbaijan 0
Bahamas 0
Bahrain 0
Bangladesh 0
Barbados 0
Belarus 98
Belgium 98
Belize 0
Benin 0
Bhutan 0
Bolivia 0
Bosnia and Herzegovina 98
Botswana 0
Brazil 0
Brunei 0
Bulgaria 98
Burkina Faso 0
Burma 0
Burundi 0
Cabo Verde 0
Cambodia 0
Cameroon 0
Canada 0
Central African Republic 0
Chad 0
Chile 0
China 0
Colombia 0
Congo (Brazzaville) 0
Congo (Kinshasa) 0
Costa Rica 0
Cote d’Ivoire 0
Croatia 98
Cuba 0
Cyprus 0
Czechia 98
Denmark 196
Diamond Princess 0
Djibouti 0
Dominica 0
Dominican Republic 0
Ecuador 0
Egypt 0
El Salvador 0
Equatorial Guinea 0
Eritrea 0
Estonia 98
Eswatini 0
Ethiopia 0
Fiji 0
Finland 98
France 98
Gabon 0
Gambia 0
Georgia 0
Germany 98
Ghana 0
Greece 98
Grenada 0
Guatemala 0
Guinea 0
Guinea-Bissau 0
Guyana 0
Haiti 0
Holy See 98
Honduras 0
Hungary 98
Iceland 98
India 0
Indonesia 0
Iran 0
Iraq 0
Ireland 98
Israel 0
Italy 98
Jamaica 0
Japan 0
Jordan 0
Kazakhstan 0
Kenya 0
Kosovo 98
Kuwait 0
Kyrgyzstan 0
Laos 0
Latvia 98
Lebanon 0
Liberia 0
Libya 0
Liechtenstein 98
Lithuania 98
Luxembourg 98
Madagascar 0
Malawi 0
Malaysia 0
Maldives 0
Mali 0
Malta 0
Mauritania 0
Mauritius 0
Mexico 0
Moldova 98
Monaco 98
Mongolia 0
Montenegro 98
Morocco 0
Mozambique 0
MS Zaandam 0
Namibia 0
Nepal 0
Netherlands 98
New Zealand 0
Nicaragua 0
Niger 0
Nigeria 0
North Macedonia 98
Norway 98
Oman 0
Pakistan 0
Panama 0
Papua New Guinea 0
Paraguay 0
Peru 0
Philippines 0
Poland 98
Portugal 98
Qatar 0
Romania 98
Russia 0
Rwanda 0
Saint Kitts and Nevis 0
Saint Lucia 0
Saint Vincent and the Grenadines 0
San Marino 98
Sao Tome and Principe 0
Saudi Arabia 0
Senegal 0
Serbia 98
Seychelles 0
Sierra Leone 0
Singapore 0
Slovakia 98
Slovenia 98
Somalia 0
South Africa 0
South Korea 0
South Sudan 0
Spain 98
Sri Lanka 0
Sudan 0
Suriname 0
Sweden 98
Switzerland 98
Syria 0
Taiwan* 0
Tanzania 0
Thailand 0
Timor-Leste 0
Togo 0
Trinidad and Tobago 0
Tunisia 0
Turkey 0
Uganda 0
Ukraine 0
United Arab Emirates 0
United Kingdom 294
Uruguay 0
US 0
Uzbekistan 0
Venezuela 0
Vietnam 0
West Bank and Gaza 0
Western Sahara 0
Yemen 0
Zambia 0
Zimbabwe 0 
nrow(frec_europa)
## [1] 185
# table crea una tabla de frecuencias, con una columna llamada Freq
table(datos_europa$Pais_Region) %>% kable()
Var1 Freq
Afghanistan 0
Albania 98
Algeria 0
Andorra 98
Angola 0
Antigua and Barbuda 0
Argentina 0
Armenia 0
Australia 0
Austria 98
Azerbaijan 0
Bahamas 0
Bahrain 0
Bangladesh 0
Barbados 0
Belarus 98
Belgium 98
Belize 0
Benin 0
Bhutan 0
Bolivia 0
Bosnia and Herzegovina 98
Botswana 0
Brazil 0
Brunei 0
Bulgaria 98
Burkina Faso 0
Burma 0
Burundi 0
Cabo Verde 0
Cambodia 0
Cameroon 0
Canada 0
Central African Republic 0
Chad 0
Chile 0
China 0
Colombia 0
Congo (Brazzaville) 0
Congo (Kinshasa) 0
Costa Rica 0
Cote d’Ivoire 0
Croatia 98
Cuba 0
Cyprus 0
Czechia 98
Denmark 196
Diamond Princess 0
Djibouti 0
Dominica 0
Dominican Republic 0
Ecuador 0
Egypt 0
El Salvador 0
Equatorial Guinea 0
Eritrea 0
Estonia 98
Eswatini 0
Ethiopia 0
Fiji 0
Finland 98
France 98
Gabon 0
Gambia 0
Georgia 0
Germany 98
Ghana 0
Greece 98
Grenada 0
Guatemala 0
Guinea 0
Guinea-Bissau 0
Guyana 0
Haiti 0
Holy See 98
Honduras 0
Hungary 98
Iceland 98
India 0
Indonesia 0
Iran 0
Iraq 0
Ireland 98
Israel 0
Italy 98
Jamaica 0
Japan 0
Jordan 0
Kazakhstan 0
Kenya 0
Kosovo 98
Kuwait 0
Kyrgyzstan 0
Laos 0
Latvia 98
Lebanon 0
Liberia 0
Libya 0
Liechtenstein 98
Lithuania 98
Luxembourg 98
Madagascar 0
Malawi 0
Malaysia 0
Maldives 0
Mali 0
Malta 0
Mauritania 0
Mauritius 0
Mexico 0
Moldova 98
Monaco 98
Mongolia 0
Montenegro 98
Morocco 0
Mozambique 0
MS Zaandam 0
Namibia 0
Nepal 0
Netherlands 98
New Zealand 0
Nicaragua 0
Niger 0
Nigeria 0
North Macedonia 98
Norway 98
Oman 0
Pakistan 0
Panama 0
Papua New Guinea 0
Paraguay 0
Peru 0
Philippines 0
Poland 98
Portugal 98
Qatar 0
Romania 98
Russia 0
Rwanda 0
Saint Kitts and Nevis 0
Saint Lucia 0
Saint Vincent and the Grenadines 0
San Marino 98
Sao Tome and Principe 0
Saudi Arabia 0
Senegal 0
Serbia 98
Seychelles 0
Sierra Leone 0
Singapore 0
Slovakia 98
Slovenia 98
Somalia 0
South Africa 0
South Korea 0
South Sudan 0
Spain 98
Sri Lanka 0
Sudan 0
Suriname 0
Sweden 98
Switzerland 98
Syria 0
Taiwan* 0
Tanzania 0
Thailand 0
Timor-Leste 0
Togo 0
Trinidad and Tobago 0
Tunisia 0
Turkey 0
Uganda 0
Ukraine 0
United Arab Emirates 0
United Kingdom 294
Uruguay 0
US 0
Uzbekistan 0
Venezuela 0
Vietnam 0
West Bank and Gaza 0
Western Sahara 0
Yemen 0
Zambia 0
Zimbabwe 0 
# los paises con frecuencia no nula serían los de Europa (???)
# P: Porqué al hacer: table(datos_europa$Pais_Region) %>% kable()
# aparecen países que no están en datos_europa, y es necesario aplicar un filter(Freq > 0) .
# R: Porque la columna de país tiene TODOS los países ya que son factores, y se guardan incluso para indicar que hay cero observaciones
# filter no deja sobre class(table), así que lo convertimos a data.frame
table(datos_europa$Pais_Region) %>% 
  as.data.frame() %>%
  filter(Freq > 0) %>%
  kable()
Var1 Freq
Albania 98
Andorra 98
Austria 98
Belarus 98
Belgium 98
Bosnia and Herzegovina 98
Bulgaria 98
Croatia 98
Czechia 98
Denmark 196
Estonia 98
Finland 98
France 98
Germany 98
Greece 98
Holy See 98
Hungary 98
Iceland 98
Ireland 98
Italy 98
Kosovo 98
Latvia 98
Liechtenstein 98
Lithuania 98
Luxembourg 98
Moldova 98
Monaco 98
Montenegro 98
Netherlands 98
North Macedonia 98
Norway 98
Poland 98
Portugal 98
Romania 98
San Marino 98
Serbia 98
Slovakia 98
Slovenia 98
Spain 98
Sweden 98
Switzerland 98
United Kingdom 294 
# Sacamos una foto de como estaba Europa el día de inicio del estado de alarma en España
datos_europa %>%
  filter(Fecha == ymd("2020-03-15")) %>%
  kable()
Provincia_Estado Pais_Region Latitud Longitud Fecha Casos_Confirmados Casos_Muertos Casos_Recuperados Casos_Enfermos
Albania 41.15330 20.16830 2020-03-15 42 1 0 41
Andorra 42.50630 1.52180 2020-03-15 1 0 1 0
Austria 47.51620 14.55010 2020-03-15 860 1 6 853
Belarus 53.70980 27.95340 2020-03-15 27 0 3 24
Belgium 50.83330 4.00000 2020-03-15 886 4 1 881
Bosnia and Herzegovina 43.91590 17.67910 2020-03-15 24 0 0 24
Bulgaria 42.73390 25.48580 2020-03-15 51 2 0 49
Croatia 45.10000 15.20000 2020-03-15 49 0 1 48
Czechia 49.81750 15.47300 2020-03-15 253 0 0 253
Faroe Islands Denmark 61.89260 -6.91180 2020-03-15 11 0 0 11
Denmark 56.26390 9.50180 2020-03-15 864 2 1 861
Estonia 58.59530 25.01360 2020-03-15 171 0 1 170
Finland 64.00000 26.00000 2020-03-15 244 0 10 234
France 46.22760 2.21370 2020-03-15 4499 91 12 4396
Germany 51.00000 9.00000 2020-03-15 5795 11 46 5738
Greece 39.07420 21.82430 2020-03-15 331 4 8 319
Holy See 41.90290 12.45340 2020-03-15 1 0 0 1
Hungary 47.16250 19.50330 2020-03-15 32 1 1 30
Iceland 64.96310 -19.02080 2020-03-15 171 5 8 158
Ireland 53.14240 -7.69210 2020-03-15 129 2 0 127
Italy 43.00000 12.00000 2020-03-15 24747 1809 2335 20603
Latvia 56.87960 24.60320 2020-03-15 30 0 1 29
Liechtenstein 47.14000 9.55000 2020-03-15 4 0 0 4
Lithuania 55.16940 23.88130 2020-03-15 12 0 1 11
Luxembourg 49.81530 6.12960 2020-03-15 59 1 0 58
Moldova 47.41160 28.36990 2020-03-15 23 0 0 23
Monaco 43.73330 7.41670 2020-03-15 2 0 0 2
Montenegro 42.50000 19.30000 2020-03-15 0 0 0 0
Netherlands 52.13260 5.29130 2020-03-15 1135 20 2 1113
North Macedonia 41.60860 21.74530 2020-03-15 14 0 1 13
Norway 60.47200 8.46890 2020-03-15 1221 3 1 1217
Poland 51.91940 19.14510 2020-03-15 119 3 0 116
Portugal 39.39990 -8.22450 2020-03-15 245 0 2 243
Romania 45.94320 24.96680 2020-03-15 131 0 9 122
San Marino 43.94240 12.45780 2020-03-15 101 5 4 92
Serbia 44.01650 21.00590 2020-03-15 48 0 0 48
Slovakia 48.66900 19.69900 2020-03-15 54 0 0 54
Slovenia 46.15120 14.99550 2020-03-15 219 1 0 218
Spain 40.00000 -4.00000 2020-03-15 7798 289 517 6992
Sweden 63.00000 16.00000 2020-03-15 1022 3 1 1018
Switzerland 46.81820 8.22750 2020-03-15 2200 14 4 2182
Channel Islands United Kingdom 49.37230 -2.36440 2020-03-15 3 0 0 3
Isle of Man United Kingdom 54.23610 -4.54810 2020-03-15 0 0 0 0
United Kingdom 55.37810 -3.43600 2020-03-15 1140 21 18 1101
Kosovo 42.60264 20.90298 2020-03-15 0 0 0 0

\[ Distancia Euclidea: d(x, y) = \sqrt{(x_{Lat}-y_{Lat})^2+(x_{Long}-y_{Long})^2}\]

# declaramos una función para calcula distancia euclidea entre dos puntos con coordenadas
distancia_grados = function(x, y){
  sqrt((x[1]-y[1])^2 + (x[2]-y[2])^2)
}

# funcion para la distancia en grados a Potsdam, con parametro x, del punto origen a calcular
distancia_grados_potsdam = function(x){
  potsdam = c(52.366856, 13.906734)
  distancia_grados(x, potsdam)
}

dist_potsdam = apply(cbind(datos_europa$Latitud, datos_europa$Longitud),
                     MARGIN = 1,
                     FUN = distancia_grados_potsdam)

datos_europa %<>%
  mutate(Dist_Potsdam = dist_potsdam)

# filtrar por las fechas del viaje, y un ciruclo alrededor de 4 grados
datos_europa %>%
  filter(between(Fecha, dmy("2-3-2020"), dmy("7-3-2020")),
         dist_potsdam < 4) %>%
  kable()
Provincia_Estado Pais_Region Latitud Longitud Fecha Casos_Confirmados Casos_Muertos Casos_Recuperados Casos_Enfermos Dist_Potsdam
Czechia 49.8175 15.473 2020-03-02 3 0 0 3 2.992057
Czechia 49.8175 15.473 2020-03-03 5 0 0 5 2.992057
Czechia 49.8175 15.473 2020-03-04 8 0 0 8 2.992057
Czechia 49.8175 15.473 2020-03-05 12 0 0 12 2.992057
Czechia 49.8175 15.473 2020-03-06 18 0 0 18 2.992057
Czechia 49.8175 15.473 2020-03-07 19 0 0 19 2.992057 
world <- ne_countries(scale = "medium", returnclass = "sf")

# Como 'Pais_Region' es un factor, el renombrado (siguente instrucción) da error, no existe ese factor level
# Añadimos el nivel de factor 'United States', los niveles serán la combinación de los niveles actuales
# 'levels(datos$Pais_Region)' más el nuevo nivel 'United States'
datos$Pais_Region = factor(datos$Pais_Region, levels = c(levels(datos$Pais_Region), "United States"))

# EEUU no aparece porque en 'datos' es 'US' y 'world' es 'United States'
datos[datos$Pais_Region=="US",]$Pais_Region = "United States"

# Vamos a tratar de cruzar los datos de nuestra tabla, para pintar los paises
# Hacemos un inner_join como en query, y la columna 'name' de 'world' contra 'Pais_Region' de 'datos'
world %>%
  inner_join(datos, by = c("name" = "Pais_Region")) %>%
  filter(Fecha == dmy("15-03-2020")) %>%
  ggplot() +
  geom_sf(color = "black", aes(fill = Casos_Confirmados)) + 
#  coord_sf(crs="+proj=laea +lat_0=50 +lon_0=10 +units=m +ellps=GRS80") + # Añadimos otra capa, una proyección para que el map no sea plano
  scale_fill_viridis_c(option = "plasma", trans = "sqrt") +
  xlab("Longitud") + ylab("Latitud") +
  ggtitle("Mapa del mundo", subtitle = "COVID-19") -> g
## Warning: Column `name`/`Pais_Region` joining character vector and factor,
## coercing into character vector
# Uso de plotly (como mejora deja hacer zoom)
ggplotly(g)
datos %>%
  filter(Fecha == dmy("30-03-2020")) %>% 
  ggplot(aes(Longitud, Latitud)) + # Hacemos una representacion con la estética (aes), de long, latit y representado con un punto geom_point()
  geom_point(aes(size = log(Casos_Confirmados+1), colour = Casos_Muertos )) + # En lugar de ver un simple punto que el punto sea proporcional y con una representacion en color según el nº de muertos
  coord_fixed() +
  theme(legend.position = "bottom") -> g
# Aplicamos en escala logaritmica a modo de ejemplo, para hacerlo proporcional, sumamos 1 para que no se ejecute log 0
# Usamos plotly
ggplotly(g)
# Umbral (threshold)
thr = 1000

datos %>%
  filter(Fecha == dmy("05-04-2020"), 
         Casos_Confirmados > thr) %>%
  mutate(Prop_Muertos = (Casos_Muertos / Casos_Confirmados) * 100, # Nueva columna de propoción de muertos 
         Ranking = dense_rank(desc(Prop_Muertos))) %>% # Añadimos un ranking en desc
  arrange(Ranking) %>% # Ordenamos la tabla por el ranking
  head(20) %>% # Mostramos solo los X primeros
  kable()
Provincia_Estado Pais_Region Latitud Longitud Fecha Casos_Confirmados Casos_Muertos Casos_Recuperados Casos_Enfermos Prop_Muertos Ranking
Italy 43.0000 12.0000 2020-04-05 128948 15887 21815 91246 12.320470 1
Algeria 28.0339 1.6596 2020-04-05 1320 152 90 1078 11.515151 2
France 46.2276 2.2137 2020-04-05 70478 8078 16183 46217 11.461733 3
United Kingdom 55.3781 -3.4360 2020-04-05 47806 4934 135 42737 10.320880 4
Netherlands 52.1326 5.2913 2020-04-05 17851 1766 250 15835 9.893003 5
Spain 40.0000 -4.0000 2020-04-05 131646 12641 38080 80925 9.602267 6
Indonesia -0.7893 113.9213 2020-04-05 2273 198 164 1911 8.710955 7
Belgium 50.8333 4.0000 2020-04-05 19691 1447 3751 14493 7.348535 8
Morocco 31.7917 -7.0926 2020-04-05 1021 70 76 875 6.856024 9
Egypt 26.0000 30.0000 2020-04-05 1173 78 247 848 6.649616 10
Iran 32.0000 53.0000 2020-04-05 58226 3603 19736 34887 6.187957 11
Sweden 63.0000 16.0000 2020-04-05 6830 401 205 6224 5.871157 12
Ecuador -1.8312 -78.1834 2020-04-05 3646 180 100 3366 4.936917 13
Hubei China 30.9756 112.2707 2020-04-05 67803 3210 63945 648 4.734304 14
Dominican Republic 18.7357 -70.1627 2020-04-05 1745 82 17 1646 4.699140 15
Philippines 13.0000 122.0000 2020-04-05 3246 152 64 3030 4.682686 16
Mexico 23.6345 -102.5528 2020-04-05 2143 94 633 1416 4.386374 17
Brazil -14.2350 -51.9253 2020-04-05 11130 486 127 10517 4.366577 18
Greece 39.0742 21.8243 2020-04-05 1735 73 78 1584 4.207493 19
Denmark 56.2639 9.5018 2020-04-05 4369 179 1327 2863 4.097047 20 
# Añadimos dos columnas
# Cortamos estos datos continuos (lat y long) en trozos o grupos con 'breaks'
# Para calcular el numero de divisiones, en este caso usalos la regla de Scott y Sturges para probar, luego lo cambiamos
datos$lat_class = cut(datos$Latitud, 
                      breaks =seq(from = -90, to = 90, by = 10))
                      #breaks = nclass.scott(datos$Latitud))
datos$long_class = cut(datos$Longitud,
                      breaks = seq(from = -180, to = 180, by = 10)) 
                      #breaks = nclass.Sturges(datos$Longitud))
tt = table(datos$lat_class, datos$long_class)
# Damos la vuelta a la tabla, porque el sur (lat y long negativas) aparecen arriba y el norte abajo
tt = tt[nrow(tt):1, ]
# Creamos un diagrama de mosaico, transponemos (t) el gráfico porque hay mas detalle en columnas que filas
# El gráfico esta mostrando los paises
mosaicplot(t(tt), shade = TRUE)

Análisis de datos temporal

# Análisis descriptivo de como han evolucionado los datos a lo largo del tiempo
# Agregamos datos en función de la fecha
# Agrega tres columnas (cbin), en función de la fecha, con el dataset (data), la función (FUN)
datos_por_fecha = aggregate(
  cbind(Casos_Confirmados, Casos_Muertos, Casos_Recuperados) ~ Fecha,
  data = datos,
  FUN = sum
)
datos_por_fecha$Casos_Enfermos = datos_por_fecha$Casos_Confirmados - 
  datos_por_fecha$Casos_Muertos - datos_por_fecha$Casos_Recuperados
head(datos_por_fecha)
##        Fecha Casos_Confirmados Casos_Muertos Casos_Recuperados Casos_Enfermos
## 1 2020-01-22               555            17                28            510
## 2 2020-01-23               654            18                30            606
## 3 2020-01-24               941            26                35            880
## 4 2020-01-25              1434            42                38           1354
## 5 2020-01-26              2118            56                51           2011
## 6 2020-01-27              2927            82                58           2787
tail(datos_por_fecha)
##         Fecha Casos_Confirmados Casos_Muertos Casos_Recuperados Casos_Enfermos
## 93 2020-04-23           2707737        190628            719350        1797759
## 94 2020-04-24           2811598        196718            769361        1845519
## 95 2020-04-25           2897619        202868            796131        1898620
## 96 2020-04-26           2972358        206568            842917        1922873
## 97 2020-04-27           3041759        211167            869418        1961174
## 98 2020-04-28           3116393        217153            902901        1996339
# Gráfico tipo barplot de los Casos_Confirmados en función de la Fecha
barplot(Casos_Confirmados ~ Fecha, data = datos_por_fecha)

plot(Casos_Confirmados ~ Fecha, data = datos_por_fecha, col = "blue", type = "l", main = "Casos_documentados por día en todo el mundo", xlab = Fecha, ylab ="Número de personas", log = "y")
lines(Casos_Muertos ~ Fecha, data = datos_por_fecha, col = "red")
lines(Casos_Recuperados ~ Fecha, data = datos_por_fecha, col = "green")

legend("topleft", c("Confirmados", "Muertos", "Recuperados"),
       col = c("blue", "red", "green"), pch = 1, lw = 2)

# Utilizando la librería 'xts'
# Crea un objeto xts, que el nombre de la fila es la Fecha y los valores: Casos_Confirmados, Muertos, Recuperados y Enfermos
datos_por_fecha_ts <- xts(x = datos_por_fecha[ ,2:5],
                          order.by = datos_por_fecha$Fecha)
# Con la libreria dygraphs asociada a xts pintamos un diagrama, con diferentes opciones
dygraph(datos_por_fecha_ts) %>%
  dyOptions(labelsUTC = TRUE, labelsKMB = TRUE,
            fillGraph = TRUE, fillAlpha = 0.05, 
            drawGrid = FALSE, colors = "red") %>%
  dyRangeSelector() %>%
  dyCrosshair(direction = "vertical") %>%
  dyHighlight(highlightCircleSize = 5, highlightSeriesBackgroundAlpha = 0.2,
              hideOnMouseOut = FALSE) %>%
  dyRoller(rollPeriod = 2)

Análisis datos España

# Filtramos Spain, pero además seleccionamos (select) las columnas que aportan algo
# En este caso, la columna Fecha, y aquellas que empiezan por "Casos_"
# La instrucción 'select' también condiciona el orden de las columnas
datos_spain = datos %>% 
  filter(Pais_Region == "Spain") %>%
  select(Fecha, starts_with("Casos_"))

plot(x = datos_spain$Fecha, y = datos_spain$Casos_Confirmados, main = "Casos Confirmados en España", type = "s", col = "blue", lwd = 2) # lwd = line width (ancho de la línea), type ="s" --> tipo escalón

# Utilizando la librería 'xts'
# Crea un objeto xts, que el nombre de la fila es la Fecha y los valores: Casos_Confirmados, Muertos, Recuperados y Enfermos
datos_spain_ts <- xts(x = datos_spain[ ,2:5],
                          order.by = datos_por_fecha$Fecha)
# Con la libreria dygraphs asociada a xts pintamos un diagrama, con diferentes opciones
dygraph(datos_spain_ts) %>%
  dyOptions(labelsUTC = TRUE, labelsKMB = TRUE,
            fillGraph = TRUE, fillAlpha = 0.05, 
            drawGrid = FALSE, colors = "red") %>%
  dyRangeSelector() %>%
  dyCrosshair(direction = "vertical") %>%
  dyHighlight(highlightCircleSize = 5, highlightSeriesBackgroundAlpha = 0.2,
              hideOnMouseOut = FALSE) %>%
  dyRoller(rollPeriod = 2)
# Gráfico barplot, tenemos que transponer (t), de la columna 3 a la 5
barplot(as.matrix(t(datos_spain[ ,3:5])),
        names = datos_spain$Fecha, 
        col = c("red", "green", "yellow"),
        main = "Estudio de casos por tipo en España",
        xlab = "Fecha", ylab = "Número de personas")
legend("topleft", c("Muertos", "Recuperados", "Enfermos"),
       col = c("red", "green", "yellow"), lw = 2, pch = 1
       )

# lag retrasa los datos, según n, en este caso retrasamos un día
# lead adelanta los datos, según n, en este caso adelantamos un día
# Es decir obtenemos el número de Casos diario, entre el acumulado de un día y el acumulado del día anterior
datos_spain %<>% 
  mutate(Nuevos_Casos_Confirmados = Casos_Confirmados - lag(Casos_Confirmados, n = 1), 
         Nuevos_Casos_Muertos = Casos_Muertos - lag(Casos_Muertos, n = 1),
         Nuevos_Casos_Recuperados = Casos_Recuperados - lag(Casos_Recuperados, n = 1)
         )

datos_spain %>% tail(20)
##         Fecha Casos_Confirmados Casos_Muertos Casos_Recuperados Casos_Enfermos
## 79 2020-04-09            153222         15447             52165          85610
## 80 2020-04-10            158273         16081             55668          86524
## 81 2020-04-11            163027         16606             59109          87312
## 82 2020-04-12            166831         17209             62391          87231
## 83 2020-04-13            170099         17756             64727          87616
## 84 2020-04-14            172541         18056             67504          86981
## 85 2020-04-15            177644         18708             70853          88083
## 86 2020-04-16            184948         19315             74797          90836
## 87 2020-04-17            190839         20002             74797          96040
## 88 2020-04-18            191726         20043             74797          96886
## 89 2020-04-19            198674         20453             77357         100864
## 90 2020-04-20            200210         20852             80587          98771
## 91 2020-04-21            204178         21282             82514         100382
## 92 2020-04-22            208389         21717             85915         100757
## 93 2020-04-23            213024         22157             89250         101617
## 94 2020-04-24            219764         22524             92355         104885
## 95 2020-04-25            223759         22902             95708         105149
## 96 2020-04-26            226629         23190            117727          85712
## 97 2020-04-27            229422         23521            120832          85069
## 98 2020-04-28            232128         23822            123903          84403
##    Nuevos_Casos_Confirmados Nuevos_Casos_Muertos Nuevos_Casos_Recuperados
## 79                     5002                  655                     4144
## 80                     5051                  634                     3503
## 81                     4754                  525                     3441
## 82                     3804                  603                     3282
## 83                     3268                  547                     2336
## 84                     2442                  300                     2777
## 85                     5103                  652                     3349
## 86                     7304                  607                     3944
## 87                     5891                  687                        0
## 88                      887                   41                        0
## 89                     6948                  410                     2560
## 90                     1536                  399                     3230
## 91                     3968                  430                     1927
## 92                     4211                  435                     3401
## 93                     4635                  440                     3335
## 94                     6740                  367                     3105
## 95                     3995                  378                     3353
## 96                     2870                  288                    22019
## 97                     2793                  331                     3105
## 98                     2706                  301                     3071
# Gráfico
#plot(datos_spain$Fecha, datos_spain$Nuevos_Casos_Confirmados, type = "l", col = "blue", xlab = "Fecha", y_lab = "Nuevos casos", main = "Nuevos registros en España")
#lines(datos_spain$Nuevos_Casos_Muertos, type = "l", col = "red")
#lines(datos_spain$Nuevos_Casos_Recuperados, type = "l", col = "green")

# Otra forma de definirlo
#plot(Nuevos_Casos_Confirmados ~ Fecha, data = datos_spain,
#     type = "l", col = "blue", 
#     xlab = "Fecha", ylab = "Nuevos casos", 
#     main = "Nuevos registros en España")
#lines(Nuevos_Casos_Muertos ~ Fecha, data = datos_spain, type = "l", col = "red")
#lines(Nuevos_Casos_Recuperados ~ Fecha, data = datos_spain, type = "l", col = "green")

# La linéa de Recuperados se sala por arriba de la gráfica con los últimos datos
# El 26/04/2020 hay 22.109 Nuevos_Casos_Recuperados
# Por lo tanto ponemos esa primera columna, como plot para que dimensione el gráfico
plot(Nuevos_Casos_Recuperados ~ Fecha, data = datos_spain,
     type = "l", col = "green", 
     xlab = "Fecha", ylab = "Nuevos casos", 
     main = "Nuevos registros en España")
lines(Nuevos_Casos_Muertos ~ Fecha, data = datos_spain, type = "l", col = "red")
lines(Nuevos_Casos_Confirmados ~ Fecha, data = datos_spain, type = "l", col = "blue")

legend("topleft", c("Confirmados", "Muertos", "Recuperados"),
       col = c("blue", "red", "green"), 
       lwd = 2, pch = 1)

Análisis por Cohortes

El análisis de cohortes es el análisis de comportamiento de un segmento determinado de usuarios que comparten una característica en común en un periodo de tiempo.

# Para poder comparar países, escalamos a la fecha del primer contagio de cada páis
# El día 0 sería justo el dia anterior -1
# Filtramos todos los registros que tengan 0 Casos_Confirmados
# group_by no agrupa sino hacemos summarise, para decirle que columna usar para totalizar
primer_contagio = datos %>%
  group_by(Pais_Region) %>%
  filter(Casos_Confirmados > 0) %>%
  summarise(Primer_Contagio = min(Fecha)-1)


# Creamos una tabla, cruzando la de datos con la tabla de primer_contagio
# En el by ponemos una o varias columnas con 'c', si se llamasen distinto habría que poner el simbolo igual, con el nombre de cada columna
# Como resultado tendremos la tabla de datos, con una última columna con 'Primer_Contagio'
# Añadimos otra columna, con los dias desde el primer contagio, restando fechas
# La resta de fechas da un objeto time, mejor convertirlo a entero, para facilitar los gráficos
# Como no todos los páises tienen el detalle de Provincia_Estado, los podemos agrupar, para hacer el análisis por Pais_Region, la métrica será que ahora los Casos_xx serán igual a la suma de los casos de lo que se está agrupando
data_first =datos %>%
  inner_join(primer_contagio, by = "Pais_Region") %>%
  mutate(Dias_Desde_PC = as.numeric(Fecha - Primer_Contagio)) %>%
  filter(Dias_Desde_PC >= 0) %>%
  group_by(Dias_Desde_PC, Pais_Region) %>%
  summarise(Casos_Confirmados = sum(Casos_Confirmados),
            Casos_Muertos = sum(Casos_Muertos),
            Casos_Recuperados = sum(Casos_Recuperados),
            Casos_Enfermos = sum(Casos_Enfermos)
            )

# Creamos el gráfico con ggplot y creando estéticas de dibujo
# Son muchos datos todos los países, así que podemos elegir una lista
data_first %>%
  filter(Pais_Region %in% c("Spain", "Italy", "China", "United States", "Germany")) %>%
  ggplot(aes(x = Dias_Desde_PC, y = Casos_Confirmados)) + 
  geom_line(aes(col = Pais_Region)) +
  xlab("Días desde el primer contagio") +
  ylab("Número de personas contagiadas") +
  ggtitle("Análisis de Cohortes") +
  theme(legend.position = "none") -> g # Quitamos la leyenda porque tendría todos los países, sino filtramos

# Lo metemos en una variable 'g', para poder usarlo con ggploty y que sea interactivo
ggplotly(g)

Modelos de regresión

\[y = f(x)\]

# Origen de la declaración de la pandemia "22/02/2020"
# Añadimos una columna que será la variable independiente 'x' con el número de días 
datos_spain$Dias = as.numeric(datos_spain$Fecha - dmy("22/01/2020"))

Regresión Lineal

\[y = ax+b, a,b\in \mathbb R\]

\[min_{a,b\in\mathbb R} \sum_{i=1}^n (y_i-(ax_i+b))^2\]

# Linear Model, la formula intentará predecir:
# El número de 'Casos_Confirmados' en función '~' de los 'Dias' desde el origen de la pandemia para el dataset 'datos_spain'
mod1 <- lm(Casos_Confirmados ~ Dias, data = datos_spain)
summary(mod1)
## 
## Call:
## lm(formula = Casos_Confirmados ~ Dias, data = datos_spain)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -62035 -36159   4409  33778  60781 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -60781.0     7924.1   -7.67 1.41e-11 ***
## Dias          2441.5      141.1   17.30  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 39520 on 96 degrees of freedom
## Multiple R-squared:  0.7571, Adjusted R-squared:  0.7546 
## F-statistic: 299.3 on 1 and 96 DF,  p-value: < 2.2e-16

La llamada a r para que saque los valores de las variables, solo se ve al ejecutar el markdown \[Casos\ Confirmados = 2441.5303827 Dias + -6.0780989\times 10^{4}\] \[y = b_1 * Dias + b_0\]

# Pintamos un gráfico con los datos de x e y
plot(datos_spain$Dias, datos_spain$Casos_Confirmados)
# Añadimos la recta de regresión
abline(mod1, col = "red")

# Pintamos distribucion de los errores ajustados al modelo
plot(mod1$residuals ~ mod1$fitted.values, xlab = "Valores Ajustados (predicción)", ylab = "Residuos del modelo (errores)")

# Este gráfico muestra que el modelo (regresion lienal) no es bueno
# No es Homocedastico

# Pasamos un test (hay diferentes)
# En la libreía 'car' hay una función qqPlot
residuos = mod1$residuals

# Sobre los residuos, siguen o no una distribución normal, con la media de los residuos, y una desviación standard de los propios residuos: residuos ~ (u, d2)
qqPlot(residuos, distribution = "norm", mean = mean(residuos), sd = sd(residuos))

## [1] 57 56
# Obtendremos una distribución de los errores, el intervalo de confianza, donde mucho datos se salen, se acumulan al inicio, y en la zona media
# Los cuantiles son puntos tomados a intervalos regulares de la función de distribución de una variable aleatoria
# Los cuartiles, que dividen a la distribución en cuatro partes (corresponden a los cuantiles 0,25; 0,50 y 0,75)
# Los quintiles, que dividen a la distribución en cinco partes (corresponden a los cuantiles 0,20; 0,40; 0,60 y 0,80);
# Los deciles, que dividen a la distribución en diez partes;
# Los percentiles, que dividen a la distribución en cien partes.

Regresión exponencial

Intentar aproximar no con una recta, sino con el logaritmo de la variable dependiente se pueda aproximar con una recta \[log(y) = ax+b, a,b \in \mathbb R\] si depejamos ‘y’ queda: \[y = e^{ax+b} = m e^{ax}\] Yo: para aclarar lo de m y el uso de exp(b) \[y = e^{ax+b} = e^b \cdot e^{ax} = m e^{ax}\]

# Creamos el seguno modelo, tenemos cuidado de no incluir 0, dado que log(0) no existe
# Cogemos solo las filas con datos > 0 y todas las columnas
# El número de 'log(Casos_Confirmados)' en función '~' de los 'Dias' desde el origen de la pandemia para el dataset 'datos_spain'
mod2 <- lm(log(Casos_Confirmados) ~ Dias, 
           data = datos_spain[datos_spain$Casos_Confirmados > 0, ])
summary(mod2)
## 
## Call:
## lm(formula = log(Casos_Confirmados) ~ Dias, data = datos_spain[datos_spain$Casos_Confirmados > 
##     0, ])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.73811 -1.13695  0.05916  1.27276  1.99353 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.586576   0.355696  -7.272 1.53e-10 ***
## Dias         0.182265   0.006006  30.348  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.431 on 86 degrees of freedom
## Multiple R-squared:  0.9146, Adjusted R-squared:  0.9136 
## F-statistic:   921 on 1 and 86 DF,  p-value: < 2.2e-16
# R2 es muy alto 0,91 (entre 0 y 1), su significatividad es muy alta
# Multiple R squared is simply a measure of Rsquared for models that have multiple predictor variables

\[Casos\ Confirmados = 0.0752773 \cdot e^{0.1822652 \cdot x}\]

plot(datos_spain$Dias, datos_spain$Casos_Confirmados)
# Exponencial de b * Exponencial de (a * x)
lines(exp(mod2$coefficients[1])*exp(mod2$coefficients[2]*datos_spain$Dias), col = "blue")

# Análisis de los residuos
plot(mod2$residuals ~ mod2$fitted.values, xlab = "Valores Ajustados (predicción)", ylab = "Residuos del modelo (errores)")

residuos = mod2$residuals
qqPlot(residuos, distribution = "norm", mean = mean(residuos), sd = sd(residuos))

## 98 34 
## 88 24

Modelo Potencial

Entre medias del lineal y el exponencial Se suele usar por si nos hemos pasado de frenada con el exponencial Parecido al exponencial, pero intentamos aproximar el log(y) en función del log(x) \[log(y) = a\cdot log(x)+b, a,b,\in\mathbb R\] Un número multiplicando un log es igual al log elevado al número que multiplica La exponencial del logaritmo se tachan, son funciones inversas Quedaria una constante ‘m’ por x^a \[y = e^{a\cdot log(x)+b} = e^b\cdot e^{log(x)^a} = m\cdot x^a \]

# El número de 'log(Casos_Confirmados)' en función '~' del 'log(Dias)' desde el origen de la pandemia para el dataset 'datos_spain'
# log(Dias) no nos preocupa, porque solo hay un día 0 y al filtrar por Casos_Confirmados > 0, no nos afecta porque el dia 0 no hay todavía Casos
mod3 <- lm(log(Casos_Confirmados) ~ log(Dias), 
           data = datos_spain[datos_spain$Casos_Confirmados > 0, ])
summary(mod3)
## 
## Call:
## lm(formula = log(Casos_Confirmados) ~ log(Dias), data = datos_spain[datos_spain$Casos_Confirmados > 
##     0, ])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.8834 -0.7721  0.2622  1.0007  4.6538 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -22.4552     1.1019  -20.38   <2e-16 ***
## log(Dias)     7.7310     0.2842   27.20   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.58 on 86 degrees of freedom
## Multiple R-squared:  0.8959, Adjusted R-squared:  0.8947 
## F-statistic: 739.8 on 1 and 86 DF,  p-value: < 2.2e-16
# R2 = 0,89

\[Casos\ Confirmados = 1.7694729\times 10^{-10} \cdot Dias^{0.1822652}\]

plot(datos_spain$Dias, datos_spain$Casos_Confirmados)
lines(exp(mod3$coefficients[1])*datos_spain$Dias^mod3$coefficients[2], col = "green")

plot(mod3$residuals ~ mod3$fitted.values,
     xlab = "Valores Ajustados", ylab="Residuos del modelo")

residuos = mod3$residuals
qqPlot(residuos, distribution = "norm", mean = mean(residuos), sd = sd(residuos))

## 11 12 
##  1  2
# CONCLUSION:
# Ninguno de los tres modelos es el bueno
# El modelo lineal es demasiado simple, la exponencial demasiado estricta crece demasiado rápido, y el potencial parece más calmada pero crece tarde
# Podríamos crear nuestro propio modelo, usando la función lm para aplicar algún tipo de transformación
# De los tres modelos, el mejor es el exponencial
# El log(Casos_Confirmados) ~ Dias da un resultado muy bueno, lo aceptamos, pero en lugar de Dias podriamos intentar añadir Dias + log(Dias), la parte del modelo exponencial + la parte del potencial 
# Añadimos más variables independientes y lo podemos complicar más con Dias^2
# La I indica que es una variable del modelo I de variable independiente y añadir más variables al modelo, hasta complicar el modelo todo lo que queramos
# Es la misma variable con transformaciones no lienales
mi_model <- lm(log(Casos_Confirmados) ~ Dias + log(Dias) + I(Dias^2) + I(Dias^3) + sqrt(Dias),
               data = datos_spain[datos_spain$Casos_Confirmados > 0, ])
summary(mi_model)
## 
## Call:
## lm(formula = log(Casos_Confirmados) ~ Dias + log(Dias) + I(Dias^2) + 
##     I(Dias^3) + sqrt(Dias), data = datos_spain[datos_spain$Casos_Confirmados > 
##     0, ])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.36401 -0.09434  0.00342  0.12260  0.70089 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.392e+02  9.389e+00   14.83   <2e-16 ***
## Dias         2.822e+01  2.182e+00   12.93   <2e-16 ***
## log(Dias)    3.012e+02  2.557e+01   11.78   <2e-16 ***
## I(Dias^2)   -1.168e-01  9.115e-03  -12.82   <2e-16 ***
## I(Dias^3)    3.247e-04  2.738e-05   11.86   <2e-16 ***
## sqrt(Dias)  -3.492e+02  2.808e+01  -12.44   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3477 on 82 degrees of freedom
## Multiple R-squared:  0.9952, Adjusted R-squared:  0.9949 
## F-statistic:  3396 on 5 and 82 DF,  p-value: < 2.2e-16
# R2 = 0,92
# Añadiendo la I -> R2 = 0,98
# Con todas las variables -> R2 = 0,99
# Comparar los modelos
# Recordad que en datos_spain, añadimos la columna Dias, con la diferencia entre la Fecha y el día de inicio de la pandemia "22/01/2020"

# Establecemos fecha de inicio y fin
start_date = ymd("2020-01-22")
end_date = ymd("2020-04-30")

# Generamos un vector con el número de día desde la fecha inicio hasta la fecha fin
# Las fechas a considerar serían las desde la de inicio hasta la final en saltos de 1 día
# Sumamos 1 al inicio para que no haya el 0, por lo de los log(0)
dates = seq(start_date + 1, end_date, by = "1 day")
# Días desde el origen de la pandemia
# A las fechas anteriores le restamos la de origen de la pandemia -> obtendremos una secuencia
days_since_start = as.numeric(dates - start_date)

# Ya podemos hacer las predicciones
# Como en los modelos solo hemos usado la columna Dias en las variables independientes, vamos a crear un data.frame donde la columna Dias viene dada por days_since_start
new_data = data.frame(Dias = days_since_start)

# Predicción 1 es con el modelo 1 lineal
pred1 = predict(mod1, newdata = new_data)
# Predicción 2 es con el modelo 2 exponencial
pred2 = exp(predict(mod2, newdata = new_data))
# Predicción 3 es con el modelo 3 potencial, también hay un log de la variable independiente
pred3 = exp(predict(mod3, newdata = new_data))
# Predicción 4 es con el modelo 2 exponencial
pred4 = exp(predict(mi_model, newdata = new_data))

# Creamos objeto de tipo xts
# Hay que rellenar con valores desconocidos, porque las predicciones tiene 100 elementos y los datos reales Casos_Confirmados en spain tenemos 98
datos_por_fecha_ts = xts(x = data.frame(Real = c(datos_spain$Casos_Confirmados, rep(NA, length(pred1) - length(datos_spain$Casos_Confirmados))),
                                        Mod_Lin = pred1,
                                        #Mod_Exp = pred2, # Quitamos el Exp porque nos destroza la gráfica y no vemos nada
                                        Mod_Pot = pred3,
                                        Mod_Mixt = pred4),
                         order.by = dates) # dates es variable calculada antes con la lista de fechas desde inicio hasta fin con saltas de 1 día
# Creamos el gráfico       
# Podemos ver que el modelo Mixto el inventado, está muy cerca del Real
dygraph(datos_por_fecha_ts)
# DUDA de lo anterior
# En el order.by -> 'dates' son fechas, y 'datos_spain' si tiene la columna Fecha y Dias (número de día desde el inicio). Pero las predicciones se han hecho solo con Dias, y no con fechas.
# ¿Cómo xts es capaz de "ordenar" los datos por fechas, y relacionar la columna de 'Casos_Confirmados', con las columnas de las predicciones 'Mod_Lin', 'Mod_Exp', etc.?
# Posible respuesta: realmente se juntan al crear el data.frame, y ahi se está la columna Dias en todos los elementos, y luego order.by solo ordena por las fechas
# Comprobación, vamos a crear aparte el data.frame, poemos [0:98], para evitar lo de añadir 'NA', y poder tambiñen añadir la columna Fecha
my_dataframe = data.frame(Real = datos_spain$Casos_Confirmados, 
                          Mod_Mixt = pred4[0:98], datos_spain$Fecha)
# Respuesta: data.frame solo junta columnas con el mismo número de filas y xts crea un objato que el nombre de la fila es el indicado en order.by, que solo añade esa columna de tiempo como índice, no trata de agrupar o relacionar los datos, solo los anexa
# Podriamos poner en order.by cualquier columna con un campo fecha y con el mismo número de filas que el data.frame y xts las uniría
my_xts_dates = xts(x = my_dataframe, order.by = dates[0:98])
my_xts_fecha = xts(x = my_dataframe, order.by = datos_spain$Fecha)
dygraph(my_xts_dates)
dygraph(my_xts_fecha)
# Ejemplo, vamos a crear de nuevo la variable dates pero con otro año
my_start_date = ymd("2021-01-22")
my_end_date = ymd("2021-05-01") # ponemos un día más porque 2020 fue bisiesto
my_dates = seq(my_start_date + 1, my_end_date, by = "1 day")
my_xts_dates = xts(x = my_dataframe, order.by = my_dates[0:98])
# Se ha creado un objeto igual pero los nombres de filas son fechas de 2021
dygraph(my_xts_dates)